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Q ■ Abstract 



Based on our studies done on two-dimensional autonomous systems, forced non-autonomous 
systems and time-delayed systems, we propose a unified methodology — that uses renormahzation 
group theory — for finding out existence of periodic solutions in a plethora of nonlinear dynamical 



^ ■ systems appearing across disciplines. The technique will be shown to have a non-trivial ability of 

l/^ , classifying the solutions into limit cycles and periodic orbits surrounding a center. Moreover, the 
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CN \ methodology has a definite advantage over linear stability analysis in analyzing centers. 
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I. INTRODUCTION 

The study of nonlinear differential equations [1, 2] in a two dimensional dynamical system 
is of considerable interest to researchers across disciplines. Various methods of obtaining 
approximate analytic solutions have been formulated over the years like Lindstedt-Poincare 
method, harmonic balance etc. More recently a number of new methods have been proposed, 
e.g. nonperturbative method [3, 4], (5-method [5], homotopy perturbation method [6], vari- 
ational iteration methods [7] etc. The use of renormalization group (RG) in the analysis of 
nonlinear dynamical problems [8] has been pioneered by Goldenfeld and co-workers [9-11]. 
The traditional perturbative theory (for example, the multiple scale method) relies on one's 
ability to recognize the correct scales. The use of RG on the direct perturbation expansion, 
eliminated the necessity of recognizing the correct scales — the scales emerged automatically 
on implementation of the RG. 

What we show in this paper is how this RG technique can be of use in devising a method- 
ology (described later in this paper) capable of distinguishing between different types of 
periodic solutions — centers and limit-cycles — in two-dimensional autonomous dynamical 
systems of the general form: x = P{x,y), y = Q{x,y). The presence of limit cycles in a 
model facilitates explanation of self-sustained oscillations. Limit cycles appear in wide vari- 
ety of modern researches in many fields like quantum physics, chemical physics, biophysics, 
material sciences, ecology etc. — see e.g., [12-16] respectively for recent examples. Conse- 
quently, finding out variety of methods [3, 17, 18] (however, see also [19]) for determining 
limit cycles in nonlinear problems is of current research interest. These orbits are essentially 
nonlinear in nature and occur isolated in a phase-space unlike the family of periodic orbits 
around a center. We shall also illustrate that this very technique can help one ascertain if 
a fixed point is a focus or a center. It is worth mentioning that distinguishing between a 
center and a focus (known as center problem) is one of the main and oldest problems in 
two-dimensional dynamical systems. Our method of distinguishing focus, center and limit- 
cycle will be easily shown to be extendable to two-dimensional non-autonomous systems 
and also to the extremely important class of time-delayed dynamical systems. Moreover, 
our technique will yield the correct nature of a fixed point of a nonlinear dynamical system 
when the linearization about the fixed point gives a completely wrong idea regarding the 
true nature of the fixed point. 



Attempts in solving an ordinary differential equation of the form: x + u'^x = eF{x,x), using 
a naive expansion, x{t) = xq + exi + e'^X2 + ■ ■ ■ , results in breakdown of the perturbation 
theory at times t such that e{t — tg) > 1 (where to is the initial time) due to the presence 
of secular terms. How does one apply the RG principle to the problem? We begin by ob- 
serving that a periodic solution can be expressed as a Fourier series with amplitude A and 
phase of the lowest harmonic, determining the amplitude and phase of the higher order 
ones. The amplitude and phase are quantities that will flow. To regularize the perturbation 
series, RG technique first introduces an arbitrary time r with a view to splitting t — to as 
{t — t) + {t — to) and absorbing the terms containing r — to into the respective renormalized 
counterparts A and of Aq and 0o. Aq and 0o are the constants of integration determined 
at to. This is completely similar to divergence in field theories where a physical quantity 
{e.g., two point correlation function) diverges as the cutoff A — )■ oo. If we are discussing a 
physical variable, then the answer has to be finite and while this is achieved in field theory 
by constructing running coupling constants, it is done for the differential equation by intro- 
ducing an arbitrary time scale r and letting the amplitude and phase depend on r. At the 
end of the process one arrives at the RG-fiow-equations for A and 0: 

So, we see that the RG naturally leads to flow equations. In this respect it is akin to 
the Bogoliubov-Krylov method [1]. But as mentioned earlier, the advantage lies in the 
fact that RG uses naive perturbation theory; and we do not need to anticipate scales (as in 
multiple scales method) or make an assumption about slowly varying amplitudes and phases 
(Bogoliubov-Krylov) . 

For autonomous systems, / and g are generally function of A alone. We propose to use flow 
equations (1.1) and (1.2) to differentiate between oscillators which are of the center variety 
and limit cycles. The center type oscillation consists of a continuous family of closed orbits 
in phase space, each orbit being determined by its own initial condition. This implies that 
the amplitude A is fixed once the initial condition is set. This must lead to 

This statement is exact and is not tied to any perturbation theory argument. For the limit 
cycle on the other hand 



and f{A) must be such that the flow has a fixed point. The fixed point has to be stable for 
the limit cycle to be stable. Also, if A = is a fixed point of equation (3), then we have a 
focus. 

This extremely simple prescription, though not proved rigorously, appeals to one's intuition 
when one notes that (i) A = means the assumed periodic solution has zero amplitude and 
hence hints at focus, (ii) f^^^ = OVA > hints at a family of non-isolated periodic orbits 
surrounding the fixed point and therefore existence of center is implied, and {in) vanishing 
of dA/dr sX A = Ai ^ Q logically indicates that an isolated periodic orbit of amplitude Ai 
happens to be surrounding the fixed point. 

The calculation of f{A) requires the use of perturbation theory. Application of perturbation 
theory is possible only if one can locate a center — this is the basic periodic state. Locating a 
center can sometimes be straightforward e.g. xi = X2, X2 = —dV/dxi, where ^ is a general 
anharmonic potential: V = xf/2 + Xix^/S + Xixf/A. Here {xi, X2) = (0, 0) is a linear center 
around which perturbation theory can be done. Similarly for the Van der Pol oscillator 
a^i = 3^2, ^2 = kxi{xl — 1) + u'^xi, the origin is a center for A; = 0. In the Lotka-Volterra 
model — ii = xi — X1X2, X2 = —X2 + a;ia;2 the origin is a saddle and (1, 1) is the center. 
Shifting the origin to the center is the first step of the process of determining the function 
f{A). In this case, of course, f{A) = since the periodic state in the predator-prey model 
is a center like state. 

A more complicated situation arises in the Belushov-Zhabotinsky reaction [20, 21] system. 
In that case, a transfer of origin to the fixed point will have to be followed by a proper 
setting of parameters in the problem to make the origin a center which is the starting point 
of all perturbation theory. This raises the problem that the given dynamical system may 
not have a relevant parameter, e.g. the well known paradigm for the limit cycle 

z = {l+i)z-/3\z\'^z (4) 

where z = x + iy is the complex variable and /3 > 0. The only fixed point is the origin and 
it is an unstable focus for all f3. We can overcome this difficulty by considering the more 
general system 

z = (ai + ia2)z — (3\z\'^z (5) 

The origin is now a stable focus for ai < 0, unstable focus for ai > and a center for 
tti = 0. It is this center about which one can set up a perturbation theory. The perturbative 



evaluation of f{A) and g{A) consequently involves the following initial steps: 

1. Find the fixed points of the system and identify linear centers. 

2. If there are no linear centers, extend the parameter space and see if a linear center can 
be located as the parameters are changed. 

3. For every linear center, thus located, we need to check the existence of a limit cycle 
by perturbatively constructing f{A) and g{A). 

II. CENTER 

In this section, we take up the study of center. The first example of the nonlinear oscillator 
will be dealt with in detail — the subsequent ones will be handled briefiy. 

A. Unforced DufRng oscillator 

The equation of motion of this damped nonlinear oscillator is 

X + kx + iJ^x + \x^ = (6) 

We notice that a linear center exists for k = X = 0. Hence the perturbation theory will have 
to be built around this limit. We expand 

X = Xq + kx[ + Xxi + k'^x'2 + \^X2 + k\x2 + . . . (7) 

Putting Eq.(7) in Eq.(6), we obtain: 

Xq + oj'^xq = (8) 

Xi + iJ^Xi = —xl (9) 

x'l + u'^x[ = —Xo (10) 

With the initial condition set as x(t = 0) = A and x{t = 0) = 0, we write the solution of 
Eq.(8) as 

Xq = AcOSUt (11) 



We note that xq picks up the initial condition and hence Xi{t = 0) = Xi{t = 0) = for all 
i > 1. Accordingly, Eq.(9) and Eq.(lO) now read: 

Xi + oj'^xi = {cos 3ujt + 3 COS ut) (12) 

x'l + oj'^x'i = uAsinut (13) 

giving rise to the following solutions respectively 

xi = —— — tsina;t+ (cos 3a;t — cos cut) (14j 

oU! 32uj 

x^ = tcosut-\ smut (15) 

2 ZLU 

keeping in mind the initial conditions. At this order, the displacement of the oscillator is 

3XA^ . AA3 ^ ^ kA kA . 
xit) = A cos ut tsinut-\ -{cos 3ujt — cos ut) 1 cos cut H smcut 

3AA3 ^ ^ AA3 

= A cos cut (t — T + t) smcut H -(cos3cut — cos cut) 

8cu 32cu^ 

— — (t — r + r) coscut + -— smcut (16) 

2 2ijj 

where we have split the interval to t as to r and r to t. To remove the divergences, we 
introduce the renormalization constants ^i(0,r) and Z2{0,t) as 

A = A{t = 0) = A{t)Zi{0, t) (17a) 

O = 0(t = O) = 0(r)+22(O,r) (17b) 

The renormalization constants have the expansion 

Zi{0,T) = l + ai\ + a[k + ... (18a) 

Z2{0,T)=bi\ + b[k + ... (18b) 

so that the constants Oj and bi can be chosen order by order to remove divergences at each 
order. In terms of A{t) and 0(t), we can write Eq.(16) as 

x{t) = A{t) [1 + aiA + a[k] cos(cut + 0(r) + biX + b[k) 

(t -t + t) sin(cut + 0) H r (cos 3(cut + 0) - cos(cut + 0)) 

8cu 32cu^ 

kA kA 
(t -t + t) cos(cut + 0) H sin(cut + 0) 

2 2uj 

= A{t) cos(cut + 0) + (aiA + a[X)A{T) cos(cut + 0) - (6iA + b[k)A{T) sin(cut + 0) 

3\A^ \A^ 

— (t -t + t) sin(cut + 0) + TT— ^(cos 3(cut + 0) 

8cu 32cu^ 

kA kA 

- cos(cut + Q)) - ^7-{t-T + T) cos(cut + 0) + 7^ sin(cut + 0) (19) 

2 2cu 



correct to C(A) and 0{k). We chose a[ = ^, ai = 0, b[ = and bi 



3A 

' 8ui 



T to write 



Eq.(20) as 



3AA3 



x{t, t) = A{t) cos{ujt + 0) — (t - r) sin(a;t + 6) + — --^(cos3(a;t + 0) - cos{oot + 0)) 



32^ 



— — (t - r) cos(cjt + 0) + 7^ sin(cjt + 0) 

2 LUJ 



(20) 



We now impose the condition that xif) has to be independent of r i.e. -— = and this 

dr 

yields (to the lowest order) 



dA kA 

rf7 ~ 2~ 
dQ 3AA2 



dr 8u 



(21a) 
(21b) 



integrating to A = Aqc '^^'^ and = 0o + ^^;;j-t- Final removal of r requires setting r = t 
and then we have 



X 



(t) = Aoe^'"'/^ 



cos 



3AA2 , 

. + ^)t + 0o 



XA 



kAn 



f — —^(cos 3(a;t + 0o) - cos{cjt + 0o)) + — -^ sin(wt + 0o 



32cj 



2a; 



(22) 



For /c = 0, we have the conservative anharmonic oscillator 



X + 00 X + Xx =0 



(23) 



for which the fixed point (0, 0) in the x — x plane [i.e. x — y plane) is a center and then as 
expected 



dA 
1^ 







(24) 



with 



XA^ 
x{t) = Aq cos nt + -— ^ [cos 3nt - COS nt] + C(A^) 



(25) 



where 



8u 



(26) 



The standard results [22] for the oscillator have, thus, been correctly captured and we find 
that the emergence of x = i; = as a center is confirmed by the fact that dA/dr = 0. It 



is worth pointing out a couple of features in the calculation. The diverging terms in the 
perturbative solution come from the lowest harmonic sine and cosine terms on the right 
hand side (inhomogeneous term). The sine term is responsible for the amplitude flow dA/dr 
and the cosine term is responsible for the phase flow. To derive the amplitude equation, 
this is what we need to concentrate on. Keeping this in mind, if we examine the structure 
of higher order terms, we find that a sine term is never generated on the right hand side 
and hence dA/dr = at all orders. This question of lowest harmonic sine and cosine 
terms is quite general and we can subsequently use this to write down the flow equation by 
inspection. It should be borne in mind that although the lowest harmonic is unnecessary 
for writing down the flow equation, it is imperative to have all the relevant harmonics for 
writing down the actual solution x{t) at any order. For the anharmonic oscillator of Eq.(23), 
the phase space trajectory is 

i^ + oo^x^ + -x^ = constant = Al + -A^ (27) 

where {Aq, 0) is the initial condition for the trajectory. It is straightforward to check that 
Eq.(25) and Eq.(26) are in exact agreement with Eq.(27) to 0{X). Thus, perturbatively the 
correct phase portrait is obtained as it should. But as is well known, the perturbation series 
is not convergent and things are bound to get worse as we go to higher values of A. This is 
a problem with all perturbative approaches. 

B. Lotka-Volterra System 

Having explained in detail the case of the anharmonic oscillator, we now turn to the 
predator-prey model which is also known to have oscillatory trajectories. The prey popula- 
tion is X and the predator population is y, with the dynamics given by 

— = X — xy (28a) 

^ = -y + xy (28b) 

The origin turns out to be a saddle and there is a center at (1, 1). We shift the origin to the 
center — a procedure which will be followed regularly in our perturbative calculation of the 



periodic trajectories. Accordingly we define 

X = X + 1 
y = Y + 1 

and write the system as: 

X = -Y -XY (29a) 

Y = X + XY (29b) 

The perturbative theory proceeds by imagining the existence of a couphng constant A in 
terms of which we have 

X = -Y - \XY (30a) 

Y = X + \XY (30b) 

We expand X and F as X = Xq + AXi + A^Xs + ■ ■ ■ and F = Fq + AFi + A^Fs + • • • 
respectively to subsequently arrive at 

Xo = -Yo (31a) 

Yo = Xo (31b) 

Xi = -Fi - XoFo (31c) 

Fi = Xi + XoFo (31d) 

X, = -Y,-{XoY, + YoX,) (31e) 

F2 = X2 + (XoFi + FoXi) (31f) 

and so on. Clearly Xq = A cost (initial condition being x = A, x = at t = 0) and 
Yq = A sin t. At the next order 

Xi+Xi = -(Lo + Lo) (32) 

where Lq = XqYq and hence 

Xi + Xi = - ( — sin2t + A2cos2t J (33) 
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FIG. 1: The initial population density is plotted against the corresponding frequency. The figure 
makes it clear the RG calculation is reasonably correct. 



There are no resonating terms on the right hand side and thus 

Xi = — (sin 2t — sin t) -\ (cos 2t — cos t) 

6 3 

keeping the initial condition Xi{t = 0) = Xi(t = 0) = in mind. 
At the succeeding order 

X2 + X2 = —{LiQ + Lio) 



(34) 



(35) 



where Liq = XqYi + YqXi. We see immediately that the resonant term on the right hand 
side of Eq.(35) is A'^ cost/12 and in keeping with our previous discussion, it follows that 



dr 

dQ _ A^ 
Ih ~ ~T2 

This gives an amplitude dependent frequency of 

A^ 



Q = 1 



12 



(36a) 
(36b) 

(37) 



which has been numerically verified as can be seen from Fig. 1. 



10 



C. A Lienard system 

A more interesting example is the dynamical system 

X = y (38a) 

y = ~x — x^ + y {1 + X + n) (38b) 

The two fixed points of the above system are at (0, 0) and (—1, 0). The former is a center and 
the latter a saddle. The oscillatory orbit needs to be investigated around (0,0). The above 
system is a second order differential equation of the Lienard variety x + xF{x) + G{x) = 0, 
where F{x) is a linear function a + f3x and G{x) = x + Xx"^. We redefine a and /3 to write 

x-kx{l + x + i2)+x + Xx'^ = (39) 

Clearly Eq.(38) is obtained for A; = 1 and A = 1. The perturbation theory has to proceed 
around the linear center which is at A; = A = 0. Accordingly, we expand 

X = xq + kxi + \x\ + k'^X2 + Ax2 + k\x2 + . . . (40) 

At different orders we have 

C>(A;°A°) : xo + xo = (41a) 

C(A;U°) : xi + xi = Xq{1 + ^i) + XqXo (41b) 

C(A;°A^) : x[+x\=xl (41c) 

0{k^\^) : X2 + a;2 = xi(l + /i) + xiio + a;oii (41d) 

C(A;°A2) : i'2 + x^ = -2xox\ (41e) 

0{k^\'^) : X2 + X2 = x[{l + /i) + aJox'i + x[xo — 2xoXi (41f) 

With the initial condition x = A, x = a.t t = 0, 

xo = A cos t (42a) 

A J^ 

xi = (1 + /i) — (t cos t - sin t) H (sin 2t - 2 sin t) (42b) 

2 6 

A? A^ J^ 

x[ = — — + -^cos2t+ — -cost (42c) 

2 6 3 

11 



The corresponding flow equations at this order are 

— = (l + /i) (43a) 

^ = (43b) 

dr 

If we were to have the possibihty of a center, then clearly /i = — 1. In this section we focus 
on the potential center and work with /i = — 1. At the next order, the flow becomes 

§ - -'4 (^^) 

— = --X^A' - -k'A' (44b) 

dr 12 24 ^ ' 

The origin is a focus if both k and A are not zero, but a center if either or both of k and A 

vanish. 

The above example is an interesting example of the usefulness of increasing the space of 

parameters in a dynamical system. Starting with Eq.(38), one would not have access to the 

different possibilities that we have been finding - e.g. the competition between the center 

and the focus. This is made possible by introducing the two parameters k and A. The 

general system has been treated in [23]. 

III. LIMIT CYCLE 

A. Van der Pol oscillator 

In this section, we consider the question of the limit cycle and we introduce the RG fiow 
by recalling the calculations of Chen et al[ll] for the Van der Pol oscillator. This system is 
represented by 

x + ex (x^ - l) + w^a; = (45) 

When looked at as the second order dynamical system x = y, y = —ty{x'^ — 1) — oo'^x, there 
is a fixed point at the origin which is a stable focus for e < and unstable focus for e > 0. 
The fixed point is a center for e = and we base the perturbation expansion around e = 0, 
expanding x as 

x{t) = xo{t) + exi{t) + e'^X2{t) + ... (46) 



12 



At different orders of e, 



Xq + OJ Xq = 

xi + u'^xi = -Xo{xl - 1] 



(47) 
(48) 



We work with initial condition x = Aq aX t = and x = at t = 0. Keeping this in mind, 
we arrive at: 



xq = Aq cos ut 

2 -^"-1 



Xi 



- { An-'-^ ] tcosiot 



32a;2 



(sinSwt — Ssinwt) 



leading to 



X = Aq cos ut + e 



- I An ] t COS ut ^ (sin 3ut — 3 sin cut) 

2 1 4 y 32a;2 ^ ' 



(49) 
(50) 

(51) 



As before we split the interval to t as to r and r to t, define the renormalization constants 
Zi and Z2 by the relation 



A = Air)Z^iO,r) 

O = 0(t = O) = 0(r)+22(O,r) 
The renormalization constants Zi and Z2 can be expanded as 

Zi{0,t) = l + aie + a2e'^ + ... 

Z2{0,T)=Pie + l32e' + ... 
To 0{e), we now have 



(52a) 
(52b) 

(53a) 
(53b) 



X 



(t) = A{1 + cue + 026^ + ...) [cos{ujt + 0) - (/3ie + /^se^ + . . . ) sin(wt + 



+ - (Ao-^Vt-r + r)cos(c^t + 0) ^^° 



4 



32w^ 



(sin 3a;t — 3 sin ut) 



We choose 



. Mn - ^ 1 r 



/3i = 



to remove divergence from the past. We are now left with 



(54) 

(55) 
(56) 



e / A^\ eA^ 

x{t) = Acos{ut + Q) + - IAq - ^ ] {t - t) cos{ut + 0) - — -^ (sin 3ut - 3 sin ut) (57) 
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We now impose the condition that dx/dr = since r is an arbitrary time and x{t) cannot 
depend on where one puts the initial condition. This leads to 

f - ° (-) 

The remaining r-dependence in x{t) is removed by setting t = t and thus 

x{t) = Acos{ujt + Q) -— -^ {sin 3ut- 3 sin ut) (60) 

The flow equation has a stable fixed point at A^ = 4 and this gives the usual Van der Pol 
limit cycle of radius 2 for small e. 



B. Lienard equation 

We now return to the Lienard equation of Eq.(39) and consider what happens at the 
second order if /x + 1 7^ 0. We find in a manner identical to that outlined above, 

dr ^^ ^ 2 8 ^ ' 

^ = (/i + 1) - -\^A^ - -k^A^ (62) 

dr ^^ ' 12 2A ^ ' 

There is a limit cycle in the system if /c 7^ and A 7^ 0; and, ^ + 1 and A have the same sign. 
The cycle is stable if k also has the same sign as /i + 1 and unstable otherwise. Opposite signs 
of k and A do not allow for the existence of a stable limit cycle — a fact easily ascertained 
by numerical experiments. The limit cycle for /c 7^ 0, A 7^ is obtained by simultaneous 
relaxing of the conditions on F{x) and G{x) in the Lienard system — the conditions being 
F{x) being either odd {center) or even {limit cycle) with G{x) odd. In our case both F{x) 
and G{x) are of mixed parity. 

C. Glycolytic oscillator 

We now turn to another example which clearly illustrates the use of shifting of origin and 
determination of the locus of Hopf bifurcation points to set up the perturbation theory and 
locate the limit cycle. This example is drawn from biology and the subject is glycolysis [25, 

14 



26]. The simplest mathematical model is that of Selkov[27] and is a 2-dimensional system. 
The variable x is the concentration of ADP (adenosine diphosphate) and y that of F6P 
(fructose-6-phosphate). The dynamics is given by 

X = —X + (a + x'^)y (63) 

y = b — {a + x^)y (64) 

where '6' is the rate of fructose production by the substrate and 'a' is the rate at which fruc- 
tose decomposes (converts to ADP). It should be noted that the presence of ADP catalyzes 
this conversion and hence 'a' is augmented io a + x^. The fixed point of the system is at 

x = h, y = b/{a + 9) (65) 

The fixed point is a stable focus for a certain parameter range and an unstable focus for 
certain others. The crossover from stable to unstable focus occurs on the boundary curve 
which is a locus of points in the a-h plane where a Hopf bifurcation occurs i.e. the fixed point 
for those values of (a, h) is a center. The curve is given by 2a = Vl + 86^ — (1 + 26^) and 
is shown in Fig. 2. For points in the shaded region the fixed point is an unstable focus and 
for these values of (a, h) a limit cycle can be shown to exist by invoking Poincare-Bendixson 
theorem. We shift the fixed point to the origin and use the new coordinates X, Y given by 

X = b + X (66) 

l/=— ^ + 1^ (67) 

To use perturbation theory, we chose (a, h) close to the boundary. Setting h = a/3/8 (the 
turning point of the curve), we take a = 1/8 — 6 to consider a point inside the boundary but 
close to it. Clearly, 6 is small and positive. To 0{6), the equation of motion reads 

X = i(X + r) + £(X,r) (68) 

y = -^X-|-£(X,F) (69) 



where 



£(X, Y) = 6{3X -Y) + J lx{X + Y)+ X^Y (70) 



We note that Eqs. (68) and (69) combine to give the oscillator. 



^ + Y = ^ (71) 



15 



n r 



T r 




FIG. 2: The curve: 2a = VT+86^ — (1 + 26^) separates the figure into shaded and unshaded 
regions. If parameters are in the shaded region, one gets limit cycle (and unstable focus) while 
unshaded region corresponds to parameters giving rise to stable focus. Linear stability analysis 
predict centers for parameters on the curve. These centers are, however, not non-linear centers. 

C has to be expanded in amplitude and the parameter 5. The amplitude will emerge as 5^^'^ 
for small 6. At the zeroth order 



Xq = A cos ( ^ + 



Yo = VScosl^ + Q + n- tan"^ v^ j 



(72) 
(73) 



The frequency is 1/v^ and the axis is tilted at an angle it — tan ^ \/2 to the X-axis. The 
amplitude A is found to from the flow which at the lowest order gives 



dr 8 

The frequency changes from the zeroth order value of —= according to the flow 



(74) 



de 

dr 



6 



A^ 



V2 ' 4V2 ''«' 

The stable fixed point A^ = 165/3 gives us the size of the limit cycle for 5 ^ 1. A typical 

small-(5 orbit is shown in Fig. 3 and bears out the correctness of the above flow. This 
technique can also be used to probe limit cycles in the more complicated model of Cera et 
al. [28]. 
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FIG. 3: Limit cycle in glycolytic oscillator for a = 0.124, b = \/0.375 and 6 = 0.001. 



D. Belushov-Zhabotinsky reaction 



An identical approach is effective for Belushov-Zhabotinsky reaction. A recent version 
[20, 21] of that reaction uses a two variable system (chlorine dioxide-iodine-malonic acid 
reaction) 



X = a — X 



y = bx [ 1 



1 + X2 

y 

l + a;2 



(76) 
(77) 



where the variable x and y are the concentrations of the intermediaries I and CIO2 which 
vary on a much faster time scale than CIO2, I2 and Malonic acid. The constants 'a' and 
'6' are parameters which depend on the rate constants and the approximately constant 
concentrations of the other reactants. We note that there is one fixed point x = a/5 and 
y = l + x'^ = l + 0^/25. Our first step is to shift the origin to (a/5, 1 + a'^/25) i.e. use the 
variables 



X 



X + 



Y + 1 + 



25 



(78) 
(79) 
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The hnear stabihty analysis of the resulting system about the fixed point X = Y = shows 
that it is a center for b = be given by 

3a 25 ,„„, 

bc = ^ 80 

5 a 

The origin is an unstable focus for b < b^. and stable for b > be- We pick a value of 'a' and 
choose b = be — S, where 6 <^ be- One carries out a perturbation analysis for the variables X 
and Y by assuming that the amplitude is small for small 6. The amplitude fiow works out 
to be 

(81) 





dA 
dr 


A*^ 4 1 


"3a^ 2 1875" 
cycle exists for positive 


where fi^ = 


(^_ 25). The limit 



[ues or 0. it IS 

apparent that as we measure the value of 'a' for which limit cycles can exist, there is a 
cyclic-fold bifurcation at a = Oc — vl9L43 — obtained by setting the expression inside 
square bracket to zero. 

E. Koch-Meinhardt reaction diffusion system 

Similar considerations apply to a model which is popular for the generation of Turing 
patterns. This is the Koch-Meinhardt reaction diffusion system [29] and for our present 
purposes only the reaction part of it is relevant. The variables x and y are the number 
densities of two species which are responsible for the pigments in the pattern and satisfy the 
reaction dynamics 

X = -x H ho- (82) 

y 

y = -y + x^ (83) 

The slowly diffusing pigment (x) is auto-catalytic and also promotes the growth of the 
antagonistic fast diffusing component (y). The rate of growth of x from the environment is 
a. The fixed point of the above system is x = 1 + a, y = {1 + a)"^. As we have made clear, 
the first step involves shifting the origin to (1 + a, (1 + cr)^), thus we define 

X = X + {l + (r) (84) 

y = Y + {l + af (85) 
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In terms of the new variables X and Y, the dynamics is 



X 



-X 



2{l + a)X-Y + X^ 



Y = 2{l + a)X -Y + X^ 



(86) 
(87) 



In perturbation theory, we are interested in the small amplitude oscillators around the origin 
and hence we can expand the denominator in Eq.(86) to write 



X 



-X 



l + a 



X - 



Y 



X' 



X 



l-a 



2(1 + a) 2(1 + (7 
Y X^ 2XY 



1 



Y 



Y' 



;i + a)2 (l + a)4 

y2 2Y^X 



+ 



l + a (l + a)2 ' 2{l + ay 

y3 



l + cx) 



l + a] 



l + a) 



{i + ay ••• 

The eigenvalues A of a linear stability analysis around the center are found from 



(88) 



A_i^l(A+l)+ 2 



l + a 



l + a 







and are seen to be 



A 



a 



± 



a^ 



(89) 



(90) 



l + a y (1 + ^)2 

For no positive a, can the origin be an unstable focus and hence it would seem there can be 
no limit cycle in the system. 

However at this point, we generalize this problem by introducing a decay constant 'a', so 
that the system (Eqs (84) and (85)) becomes 

„2 



The fixed point is now aX a + a, y 
at the fixed point, we have 

X = 



X 



[a + a) 



a 



X 

-X -\ \- a 

y 

-ay + x^ 



(91) 
(92) 



and using a coordinate system (A', 3^) centered 



-X 



2ia + a)X-ay + X^' 



[a + a] 



y 



[a + a] 



y = 2ia + a)X-ay + X^ 



(93) 
(94) 
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FIG. A: a — (J parameter space. The shaded regions are where hmit cycle solutions can occur. 

Linear stability analysis in this case about the fixed point (0, 0) shows that it is a center for 

l-a 



..^-. 



(95) 



This curve is shown in Fig. 4. The interior region has the unstable focus and hence in this 
part of the parameter space there is a limit cycle. The range of a is limited to < o" < 3— 2v^ 
and 'a' lies between and 1. Clearly the limit cycle at a = 1 and o" 7^ is ruled out as 
expected. 

F. Summarizing... 

The above example serves the purpose of establishing our main point that the existence 
of a limit cycle would mean a fiow equation of the form written down in Eq(l) and providing 
a method (although perturbative) of calculating the function f(A). In the process, we note 
the following facts: 

• If a stable limit cycle exists, then there must exist an unstable focus. 

• If there are not enough parameters in the system to tune the focus to a center, then 
the linear terms in the system can be supplied with variable coefficients which may be 
tuned to yield a center. Perturbation theory can be carried out around this center. 
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It may very well be that there exists a family of limit cycles surrounding a focus. Such a 
case is naturally taken care of in our methodology, as the RG flow equation for amplitude — 
dA/dr = f{A) — would result in more than one fixed points: If there actually are A^ limit 
cycles, then f{A) = will have A^ positive real roots. 

IV. NON-AUTONOMOUS SYSTEMS 

Having clearly illustrated how to distinguish between focus, center and limit cycle in 
a two-dimensional autonomous dynamical system, we now examine how our methodology 
fares in somewhat more complicated cases of non-autonomous systems and systems with 
time-delay. 

A. A damped driven ocsillator 

We begin with a damped driven oscillator 

X + uP'x + kx = F cos Vlt (96) 

which we write as, 

x + n'^x = -kx + F cos nt + {n'^ -uj^)x. (97) 

We treat k, F and Q^ — uj^ as small to perturb about a center (/c = F = ^2^ — w^ = ). 
Accordingly, proceeding as explained earlier, to the first order in all these small parameters, 
we obtain: 

dA kA FsinG rfO FcosG 

— = : — = hAcu (98) 

dT 2 2fi ' dr 2^A ^ ' 

where Acj = w — ^2. Since, ^2 is maintained externally, it cannot change, implying dQ/dr = 0. 
Also, existence of fixed point requires dA/dr = 0. Therefore, the fixed point corresponds to 
the amplitude A = F/[A;^-h4(Aw)^]^/^ and the phase 6 = tan~^[-A;/2(Aw)]. This is exactly 
in accordance with the literature of forced oscillators. The stable non-zero fixed point in the 
evolution of A corresponds to a limit cycle in accordance with what we have claimed has to 
happen. 
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B. Time delay equations 

Now let us consider oscillators with time delay. A linear form of such an oscillator satisfies 
the differential equation: 

x{t) + uj'^xit) + exit -td) = (99) 

Here we treat e as small and consider the perturbation of the center [e = 0). The RG flow 
equations upto 0{e) are found to be: 

dA sin utd dQ cosuta ,_^, 

— — = e ; — — = e (100) 

dr 2uj dr 2uj 

For a periodic orbit (center), our claim is: dA/dr should be trivially zero; and this yields 
utd = TT. This may be compared with the exact result that equation (99) exhibits oscillatory 
solution Aexp{it\/uj'^ — e) + c.c. when t^ = tx j ^{uj'^ — e). The frequency of the periodic 
orbit is seen to be w — ejloj + 0[e^^ in accordance with the expansion of the exact answer 
of equation (99). 

Similarly, one can study limit cycles too in weakly nonlinear time delayed systems with 
success [24]. The system we study next to illustrate that RG can be successfully implemented 
in such systems, is given by 

—^ + ax{t) + fix{t -td) = \ {x{t) - x^{t)) (101) 

where a, /3 A and are constants, A being small. The LHS of Eq.(lOl) constitutes the 
unperturbed system and the nonlinear terms in RHS will be treated as the perturbation. 
We proceed as usual with a naive expansion of the form x{f) = xq (t) + \xi (t) + \'^X2 (t) + ■ ■ ■ . 
At zeroth order, we have 

^^ + axoit) + /3xo{t -td) = (102) 

It is easy to see analytically that the above equation has an oscillatory solution given the 
following condition is satisfied. 

U = '-^^-IMB (103) 

Restricting ourselves to cases where the above condition holds we find 

Xo{t) = Aocosut (104) 
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where to = \/ f? — cP' ; (3 > a. Further we find: /3 sin cur = to and /3 cos cur = —a. At order 
C(A^), using Eq. (104) we have 

dxiit) 



dt 



+ axi{t) + /3xi{t — td) = XQ{t) — XQ{t) 



A \ n 

An COS Ut — COS 3ut 

4 / 4 



(105) 



A httle bit of algebra yields the solution for xi{t) which reads 

1 + atd , , , iotd 



Xiit) = o 1 cos cut 

To C(A^), thus, we have 

1 + C(td 

x[t) = Ao cos cut + A- — ; :^tcoscut + A- 



(1 + atdY + (cutd)=^ 



Utd 



-t sin cut 



(106) 



-t sin cut (107) 



(1 + atdY + (utdY (1 + atdY + (cutdY 

From this point onwards proceeding as described earlier we arrive at the RG fiow equation 
right upto 0{X^), given by 



dA 



XA{1 + atd) 



dr (1 + atdY + i^UY 



1--A' 



dQ 



XA{uta 



1--A' 



(lOJ 



4 J' cir (i + atdY + icutdYl ^' 
In accordance to our classification scheme we can immediately conclude that the system 
given by Eq.(lOl) exhibits limit cycle oscillations. Amplitude of the limit cycle is given by 
the stable fixed point A^ = 4/3. 



V. ADVANTAGE OVER LINEAR STABILITY ANALYSIS 



Before we conclude, let us witness how useful this RG technique is when one deals with the 
subtle cases of centers in nonlinear dynamical systems. It is very well known that linearized 
version of a nonlinear dynamical system may not reproduce qualitatively correct picture of 
the phase portrait near a fixed point. We now showcase the fact that while linearization 
of a certain nonlinear dynamical system wrongly establishes a fixed point as center (which 
originally is a spiral node), our methodology gives correct result. Consider the following 
dynamical system: 

(109a) 



X = —y + eax{x^ + y'^) 



y = +x + eay{x^ + y"^) 



(109b) 
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Here, £ is a small positive parameter that facilitates a trial perturbative solution of the 
form: x{t) = xq + exi + 5^X2 + ■ ■ ■ • Linear stability analysis would show that the fixed 
point (0, 0) is a center for all a. It can however be easily shown [2] by making use of polar 
coordinates, in system (109), the origin is a stable spiral when a < and an unstable spiral 
for positive a. Now, applying the RG methodology prescribed in this paper, one arrives at 
the following flow equations, upto 0{e'^): 

dA/dr = aA^; dB/dr = (110) 

One immediately notes that in accordance with our scheme of classifying focus and center, 
from the above flow equations, one can easily extract the correct information regarding the 
nature of the flxed point in system (109): if a = 0, dA/dr = VA implying that the origin 
is a center; whereas if a 7^ 0, dA/dr = iff A = 0, making the origin a focus. 
One may recall that the flxed point given by expression (65) for the glycolytic oscillator 
deflned by equations (63) and (64) was found to be a center on curve: 2a = vT+~86^— (1 + 
26^). However, this is a result of linear stability analysis where, by dint of very nature of the 
analysis technique, the flxed point is shielded from full bombardment of non-linear terms. 



For a speciflc value of (a, b) = (1/8, y 3/8) lying on the curve, one obtains the flow equation 
[setting 6 = in relation (74)]: dA/dr = — 3A^/8 (to the lowest order). One shouldn't be 
confused to observe that in accordance with the form of this flow equation, our prescription 
claims (a, 6) = (1/8, a/3/8) is actually a focus and not a center. Numerical simulations 
easily confirm this fact. 

VI. CONCLUSIONS 

To conclude, we again emphasize that this paper introduces a simple yet powerful method- 
ology — based on perturbative renormalization group theory — of identifying and classifying 
a periodic solution (limit-cycle or orbit around center) in various types of two-dimensional 
nonlinear dynamical system. This very technique can also distinguish between a focus and 
a center. The different types of two-dimesional systems that can be handled using this 
methodology include not only simpler autonomous systems but also forced non-autonomous 
systems and time-delayed systems.. Also, it has been shown that our technique yields the 
correct nature of the fixed point of a nonlinear dynamical system when the linearization 
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about it gives a completely wrong idea regarding its true nature. 

Given the inter-disciplinary nature of the subject of nonlinear dynamics and the wide re- 
search interest in investigating periodic solutions, our method should be of direct interest 
and practical use to researchers across scientific disciplines. 
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